Friend Module modProbability

    ''' <summary>
    ''' Calculate some basic statistics from a dataset
    ''' </summary>
    ''' <param name="Values">Array of numbers</param>
    ''' <param name="Sum">Sum of all values</param>
    ''' <param name="Mean">Computed mean value</param>
    ''' <param name="StdDev">Computed standard deviation</param>
    ''' <param name="Count">Number of values</param>
    ''' <param name="MinValue">Minimum value</param>
    ''' <param name="MaxValue">Maximum value</param>
    Friend Sub CalcStats(ByVal Values() As Double, ByRef Sum As Double, ByRef Mean As Double, ByRef StdDev As Double, ByRef Count As Integer, ByRef MinValue As Double, ByRef MaxValue As Double)
        Try
            Sum = 0
            Dim sum2 As Double = 0
            MinValue = Double.MaxValue
            MaxValue = Double.MinValue
            For i As Integer = 0 To Values.Length - 1
                Sum += Values(i)
                sum2 += Values(i) ^ 2
                MinValue = Math.Min(MinValue, Values(i))
                MaxValue = Math.Max(MaxValue, Values(i))
            Next

            If Values.Length = 0 Then
                Mean = 0
                StdDev = 0
                MinValue = 0
                MaxValue = 0
            ElseIf Values.Length = 1 Then
                Mean = Values(0)
                StdDev = 0
            Else
                Mean = Sum / Values.Length
                StdDev = System.Math.Sqrt((Values.Length * sum2 - Sum ^ 2) / (CDbl(Values.Length) * (CDbl(Values.Length) - 1)))
            End If
            Count = Values.Length
        Catch ex As Exception
            ErrorMsg(, ex)
        End Try
    End Sub

    ''' <summary>
    ''' Calculate stats and lower and upper fence values (used to determine outliers)
    ''' </summary>
    ''' <param name="Values">Array of numbers</param>
    Friend Sub CalcStats(ByVal Values() As Double, ByRef Mean As Double, ByRef StdDev As Double, ByRef LowerFence As Double, ByRef UpperFence As Double)
        Try
            Dim sum As Double = 0, sum2 As Double = 0
            For i As Integer = 0 To Values.Length - 1
                sum += Values(i)
                sum2 += Values(i) ^ 2
            Next

            If Values.Length <= 1 Then
                LowerFence = Double.MinValue
                UpperFence = Double.MaxValue
                Exit Sub
            End If

            Mean = sum / Values.Length
            StdDev = System.Math.Sqrt((Values.Length * sum2 - sum ^ 2) / (CDbl(Values.Length) * (CDbl(Values.Length) - 1.0)))

            Dim Q25 As Double = NormalInv(0.25) * StdDev + Mean
            Dim Q75 As Double = NormalInv(0.75) * StdDev + Mean
            Dim IQR As Double = Q75 - Q25
            LowerFence = Q25 - 1.5 * IQR
            UpperFence = Q75 + 1.5 * IQR
        Catch ex As Exception
            ErrorMsg(, ex)
        End Try
    End Sub

    ''' <summary>
    ''' Calculate percentiles in 5 percent increments and store in 19 element array (0th element is 5%); if invalid percentile is set to 0.0
    ''' </summary>
    ''' <param name="Values">Array of numbers</param>
    ''' <param name="Percentiles">19 element array of percentile values (0th element is 5%, 18th is 95%)</param>
    Friend Sub CalcStats(Values() As Double, ByRef Percentiles() As Double)
        Dim lstVals As New Generic.List(Of Double)
        lstVals.AddRange(Values)
        lstVals.Sort()
        'do percentiles at 5 percent levels, interpolating if necessary
        ReDim Percentiles(18)
        For j As Integer = 1 To 19
            Dim recNumVal As Double = j * 5.0 * (lstVals.Count + 1) / 100.0 'taken from CRC Standard Math tables
            If recNumVal < 1 Or recNumVal > lstVals.Count Then 'too few values in list; set stat to -1 (don't extrapolate)
                Percentiles(j - 1) = 0
            Else
                Dim recNum As Integer = Math.Truncate(recNumVal) 'next smaller number
                Dim result As Double = lstVals.Item(recNum - 1)
                If recNum = lstVals.Count Then 'set to highest value
                    Percentiles(j - 1) = result
                Else 'interpolate between current and next higher value
                    Percentiles(j - 1) = result + (recNumVal - recNum) * (lstVals.Item(recNum) - result)
                End If
            End If
        Next
    End Sub

    Friend Enum enumCurveFit
        Linear
        Power
        Exponential
        OneToOne
    End Enum

    Friend Structure structFitStats
        Dim XMin, XMax, YMin, YMax As Double
        Dim NumPoints As Integer
        Dim Slope As Double
        Dim Intercept As Double
        Dim CorrCoef As Double
        Dim MeanAbsError As Double
        Dim RMSError As Double
        Dim NormRMSError As Double
        Dim IndexOfAgreement As Double
    End Structure

    '#If Project = "WRDBGraph" Then

    ''' <summary>
    ''' Calculate least squares fit statistics
    ''' </summary>
    ''' <param name="Values">Two dimensional array of paired data (0 to 1, 0 to NumPts-1); X values are 0th, Y are 1st</param>
    ''' <param name="FitType">Type of least-squares fit</param>
    ''' <param name="Stats">Computed statistics</param>
    ''' <remarks>For Index of Agreement stat, assumes that x=simulated and y=observed</remarks>
    Public Function ComputeLeastSquares(Values(,) As Double, ByVal FitType As enumCurveFit, ByRef Stats As structFitStats) As Boolean
        Try
            With Stats
                .XMin = Double.MaxValue
                .XMax = Double.MinValue
                .YMin = Double.MaxValue
                .YMax = Double.MinValue

                Dim SumX As Double = 0, SumX2 As Double = 0, SumY As Double = 0, SumY2 As Double = 0, SumXY As Double = 0
                Dim SumError As Double = 0, SumError2 As Double = 0, SumIADenom As Double = 0

                .NumPoints = Values.GetLength(1)
                For i As Integer = 0 To .NumPoints - 1
                    Dim X As Double = Values(0, i)
                    Dim Y As Double = Values(1, i)
                    .XMin = Math.Min(.XMin, X)
                    .XMax = Math.Max(.XMax, X)
                    .YMin = Math.Min(.YMin, Y)
                    .YMax = Math.Max(.YMax, Y)

                    Select Case FitType
                        Case enumCurveFit.Linear, enumCurveFit.OneToOne
                        Case enumCurveFit.Power
                            If X <= 0 Or Y <= 0 Then ErrorMsg("Cannot fit logarithmic data to non-positive data") : Return False
                            X = Math.Log(X)
                            Y = Math.Log(Y)
                        Case enumCurveFit.Exponential
                            If Y <= 0 Then ErrorMsg("Cannot fit logarithmic data to non-positive data") : Return False
                            Y = Math.Log(Y)
                    End Select
                    SumX += X
                    SumX2 += X ^ 2
                    SumY += Y
                    SumY2 += Y ^ 2
                    SumXY += X * Y
                    SumError += Math.Abs(Y - X)
                    SumError2 += (Y - X) ^ 2
                    Application.DoEvents()
                Next i

                .MeanAbsError = SumError / .NumPoints
                .RMSError = Math.Sqrt(SumError2 / .NumPoints)
                .NormRMSError = .RMSError / Math.Sqrt(SumXY / .NumPoints)

                If .NumPoints > 1 And .YMax - .YMin > 0 And .XMax - .XMin > 0 Then
                    .Slope = (SumXY - SumX * SumY / .NumPoints) / (SumX2 - SumX ^ 2 / .NumPoints)
                    .Intercept = (SumY - .Slope * SumX) / .NumPoints
                    Dim XBar As Double = SumX / .NumPoints
                    Dim YBar As Double = SumY / .NumPoints
                    SumX2 = 0
                    SumY2 = 0
                    SumXY = 0
                    For i As Integer = 0 To .NumPoints - 1
                        Dim X As Double = Values(0, i)
                        Dim Y As Double = Values(1, i)
                        Select Case FitType
                            Case enumCurveFit.Linear, enumCurveFit.OneToOne
                            Case enumCurveFit.Power
                                X = Math.Log(X)
                                Y = Math.Log(Y)
                            Case enumCurveFit.Exponential
                                Y = Math.Log(Y)
                        End Select
                        SumX2 += (X - XBar) ^ 2
                        SumY2 += (Y - YBar) ^ 2
                        SumXY += (X - XBar) * (Y - YBar)
                        SumIADenom += (Math.Abs(X - YBar) + Math.Abs(Y - YBar)) ^ 2 'assumes that X is simulated; Y is measured
                    Next i
                    .CorrCoef = Math.Abs(SumXY / System.Math.Sqrt(SumX2 * SumY2))
                    .IndexOfAgreement = 1.0 - (SumError2 / SumIADenom)
                    Return True
                Else
                    Return False
                End If
            End With
        Catch ex As Exception
            ErrorMsg(, ex)
            Return False
        End Try
    End Function

    '#End If

    ''' <summary>
    ''' Compute Log Pearson Type III Frequency Factor
    ''' </summary>
    ''' <param name="Prob">Cumulative probability (decimal)</param>
    ''' <param name="Skew">Skewness coefficient</param>
    ''' <returns>Frequency factor (K)</returns>
    ''' <remarks>
    ''' Note: approximation formula taken from Handbook of Hydrology by Maidment
    ''' the approximation is accurate to ~1% at skews up to about +-2.0
    ''' At larger skews, accuracy sharply drops
    '''</remarks>
    Friend Function LP3(ByVal Prob As Double, ByVal Skew As Double) As Double
        Try
            Dim z As Double = NormalInv(Prob)
            If Skew = 0 Then
                Return z
            Else
                Return 2 / Skew * (1 + Skew * z / 6 - Skew ^ 2 / 36) ^ 3 - 2 / Skew
            End If
        Catch ex As Exception
            ErrorMsg(, ex)
        End Try
    End Function

    ''' <summary>
    ''' Standard normal distribution routine
    ''' </summary>
    ''' <param name="z">Standard normal deviate</param>
    ''' <returns>Cumulative probability (decimal)</returns>
    ''' <remarks>Requires function ERF</remarks>
    Friend Function Normal(ByVal z As Double) As Double
        Try
            If z < 0 Then
                Return 1 - 0.5 * (1.0# + ERF(System.Math.Abs(z) * System.Math.Sqrt(0.5)))
            Else
                Return 0.5 * (1.0# + ERF(System.Math.Abs(z) * System.Math.Sqrt(0.5)))
            End If
        Catch ex As Exception
            ErrorMsg(, ex)
        End Try

    End Function

    ''' <summary>
    ''' Inverse normal distribution routine
    ''' </summary>
    ''' <param name="Prob">Probability</param>
    ''' <returns>Standard Normal Deviate, Z</returns>
    ''' <remarks>Note: approximation formula taken from Handbook of Hydrology by Maidment</remarks>
    Friend Function NormalInv(ByVal Prob As Double) As Double
        Try
            If Prob <= 0 Or Prob >= 1 Then
                'MsgBox "Illegal value for probabability passed to NormalInv function: " & Prob, , "Error"
                '  Err.Raise 1050, , "Illegal value for probabability passed to NormalInv function: " & Prob
                Exit Function
            End If

            If Prob = 0.5 Then
                Return 0
            Else
                Dim y As Double
                If Prob < 0.5 Then
                    y = -System.Math.Log(2 * Prob)
                Else
                    y = -System.Math.Log(2 * (1 - Prob))
                End If
                Dim z As Double = -System.Math.Sqrt(y * y * ((4 * y + 100) * y + 205) / (((2 * y + 56) * y + 192) * y + 131))
                If Prob > 0.5 Then z = -z
                Return z
            End If
        Catch ex As Exception
            ErrorMsg(, ex)
        End Try
    End Function

    ''' <summary>
    ''' Error function
    ''' </summary>
    ''' <param name="X">Function argument</param>
    ''' <returns>ERF function value</returns>
    Function ERF(ByVal X As Double) As Double
        Dim Eps, X2 As Double
        Dim XErr, Er, r, C0 As Double
        Dim k As Short

        Eps = 0.00000000001

        X2 = X * X

        If (System.Math.Abs(X) < 3.5) Then
            Er = 1.0#
            r = 1.0#

            For k = 1 To 50
                r = r * X2 / (k + 0.5)
                Er = Er + r
                If (System.Math.Abs(r) <= System.Math.Abs(Er) * Eps) Then Exit For
            Next k

            C0 = 2.0# / System.Math.Sqrt(Math.PI) * X * System.Math.Exp(-X2)
            XErr = C0 * Er
        Else
            Er = 1.0#
            r = 1.0#

            For k = 1 To 12
                r = -r * (k - 0.5) / X2
                Er = Er + r
            Next k

            C0 = System.Math.Exp(-X2) / (System.Math.Abs(X) * System.Math.Sqrt(Math.PI))

            XErr = 1.0# - C0 * Er

            If (X < 0.0#) Then XErr = -XErr
        End If

        Return XErr
    End Function

    ''' <summary>
    ''' Return cumulative normal distribution given standard normal deviate
    ''' </summary>
    ''' <param name="Z">Standard normal deviate</param>
    ''' <returns>Cumulative probability (decimal)</returns>
    Public Function Normal_Cumulative(ByVal Z As Double) As Double
        Try
            Dim t As Double = 1 / (1 + 0.2316419 * System.Math.Abs(Z))
            t = Normal_density(System.Math.Abs(Z)) * (0.31938153 * t - 0.356563781 * t ^ 2 + 1.781477937 * t ^ 3 - 1.821255978 * t ^ 4 + 1.330274429 * t ^ 5)
            If Z > 0 Then
                Normal_Cumulative = 1 - t
            Else
                Normal_Cumulative = t
            End If
        Catch ex As Exception
            ErrorMsg(, ex)
        End Try
    End Function

    ''' <summary>
    ''' Return standard normal density value given standard normal deviate
    ''' </summary>
    ''' <param name="Z">Standard normal deviate</param>
    ''' <returns>Standard normal density value</returns>
    Public Function Normal_density(ByVal Z As Double) As Double
        Try
            Normal_density = System.Math.Exp(-Z * Z / 2) / System.Math.Sqrt(2 * Math.PI)
        Catch ex As Exception
            ErrorMsg(, ex)
        End Try
    End Function

End Module

''' <summary>
''' This module was converted from computational library purchased by Wilson Engineering
''' </summary>
Friend Module modNLFit

    Friend Function ComputeSineCurve(ByRef Date_Time() As Double, ByVal Value() As Double, ByRef CM As Double, ByRef CA As Double, ByRef CC As Double, ByRef Period As Double, ByRef StdDev As Double, ByRef Description As String) As Boolean
        Try
            Dim X() As Double
            Dim Y() As Double 'time and value
            Const NumUnknowns As Short = 3 'number of unknowns
            Dim PInit(NumUnknowns) As Double
            Dim PUpper(NumUnknowns) As Double
            Dim PLower(NumUnknowns) As Double 'following are initial guesses, upper and lower limits
            Dim PFitted(NumUnknowns) As Double
            Dim SDev(NumUnknowns) As Double 'fitted parameters and standard deviations
            Dim NumPoints, NumIter As Integer 'number of data points, number of iterations to solve
            Dim ymin, ymax As Double
            Dim Sum, dy, sum2 As Double

            'assemble arrays of data points (NLFit routine assumes base-1 arrays)

            NumPoints = 0
            ymin = Double.MaxValue
            ymax = Double.MinValue
            ReDim X(Date_Time.Length)
            ReDim Y(Value.Length)
            NumPoints = Date_Time.Length
            For i As Integer = 1 To NumPoints
                X(i) = Date_Time(i - 1)
                Y(i) = Value(i - 1)
                ymin = Math.Min(ymin, Y(i))
                ymax = Math.Max(ymax, Y(i))
            Next i

            'set up solution parameters
            PInit(1) = (ymin + ymax) / 2 : PLower(1) = ymin : PUpper(1) = ymax 'M parameter (mean value)
            PInit(2) = (ymax - ymin) / 2 : PLower(2) = 0 : PUpper(2) = (ymax - ymin) / 2 'A parameter (amplitude)
            PInit(3) = Math.PI : PLower(3) = 0
            PUpper(3) = 2 * Math.PI 'C parameter (phase angle in radians)

            'set global variable used in USERFUNC routine
            SineCurvePeriod = Period

            'call non-linear fit solution routine
            Call NLFIT(NumUnknowns, NumPoints, X, Y, PInit, PUpper, PLower, PFitted, SDev, NumIter)

            Dim Suggestions As String = "\n\nTry the following approaches to achieve an acceptable curve fit:\n\n\t\b Adjust or delete outlier data points\n\t\b Try a different curve fit (fewer parameters are easier to fit)\n\t\b Adjust your seed, minimum and maximum values.\n\t\b Try to get a good starting 'eye-ball' fit before recalculating."
            If NumIter = modNLFit.MaxIterations Then MessageBox.Show(StringFormat("The non-linear curve-fit algorithm was unable to converge on a solution in the maximum number of iterations ({0}).{1}", modNLFit.MaxIterations, Suggestions), "Non-Linear Curve Fit Non-Convergence Error", MessageBoxButtons.OK, MessageBoxIcon.Warning) : Return False

            Dim HitBoundary As Boolean
            For i As Integer = 1 To NumUnknowns
                If Math.Abs(PFitted(i) - PLower(i)) < 0.0001 Or Math.Abs(PFitted(i) - PUpper(i)) < 0.0001 Then HitBoundary = True
            Next
            If HitBoundary Or NumIter = 1 Then MessageBox.Show(StringFormat("The non-linear curve-fit algorithm stopped after {0} iterations with one or more parameters set to the minimum or maximum allowable value, indicating that an acceptable solution could not be found.{1}", NumIter, Suggestions), "Non-Linear Curve Fit Non-Convergence Error", MessageBoxButtons.OK, MessageBoxIcon.Warning) : Return False

            'compute standard deviations of the residuals
            Sum = 0
            sum2 = 0
            For i As Integer = 1 To NumPoints
                dy = (PFitted(1) + PFitted(2) * (System.Math.Sin(2 * Math.PI / Period * X(i) + PFitted(3)))) - Y(i)
                Sum = Sum + dy
                sum2 = sum2 + dy * dy
            Next i

            StdDev = System.Math.Sqrt((NumPoints * sum2 - Sum * Sum) / (CDbl(NumPoints) * (CDbl(NumPoints) - 1.0)))

            Description = String.Format("Iterations: {0}\nNo. of points: {1}\nPeriod: {2}\nMean = {3:N4}\nAmplitude = {4:N4}\nPhase Angle = {5:N4}\nStd Dev Mean = {6:N4}\nStd Dev Amp. = {7:N4}\nStd Dev Phase = {8:N4}\nStd Dev Residuals = {9:N4}", NumIter, NumPoints, Period, PFitted(1), PFitted(2), PFitted(3), SDev(1), SDev(2), SDev(3), StdDev)
            Dim s As String = String.Format("Form of the equation is:\n\n\tY = M + A * sin(2*Pi/{0} * T + C)\n\nwhere:\tM = mean value\n\tA = amplitude\n\tC = phase angle\n\tY = predicted value\n\tT = Julian date (1/1/1900 00:00 = 1.00).\n\nResults of non-linear solution fit:\n\n{1}", Period, Description)
            MessageBox.Show(StringFormat(s), "Fit Sine Curve Results", MessageBoxButtons.OK, MessageBoxIcon.Information)
            CM = PFitted(1)
            CA = PFitted(2)
            CC = PFitted(3)
            Return True
        Catch ex As Exception
            ErrorMsg(, ex)
            Return False
        End Try
    End Function


    Friend MaxIterations As Integer = 100

    ''' <summary>
    ''' This must be set before calling NLFit routine
    ''' </summary>
    Friend SineCurvePeriod As Integer

    Private Sub USERFUNC(ByRef k As Integer, ByRef X() As Double, ByRef p() As Double, ByRef Y As Double)

        ' Non-linear fitting function subroutine

        ' Input

        '  K   = index of X data array
        '  X() = array of X data points
        '  P() = array of fitting parameters

        ' Output

        '  Y = function value at X(K)

        'Fit annual sine curve to time series data (mainly temperature)

        'y = M + A [sin(bt + c)]

        'b is known (assume annual, so equals 2pi/365)
        'there are three unknowns: M, A, and C

        Static c, a, b, M As Double
        M = p(1) : a = p(2) : c = p(3)
        b = 2 * Math.PI / SineCurvePeriod 'global variable set before call
        Y = M + a * (System.Math.Sin(b * X(k) + c))
    End Sub


    'following routines taken from BNALib library; source code was purchased for WRDB project

    ' Copyright (c) 1997 C. David Eagle Jr.

    Sub DECOMPOSE(ByRef N As Integer, ByRef EMATRIX() As Double, ByRef BMATRIX(,) As Double)
        Dim Flip As Double
        Dim STETA As Double
        Dim CTETA As Double
        Dim TETA As Double
        Dim WW As Double
        Dim AC As Double
        Dim XLE As Double
        Dim WIH As Double
        Dim WIJ As Double
        Dim i As Integer
        Dim XHLONG As Double
        Dim XJLONG As Double
        Dim SCALPROD As Double
        Dim h As Double
        Dim j As Integer
        Dim KI As Integer
        Dim NI As Integer
        Dim IT As Integer
        Dim ITMAX As Integer
        Dim Eps As Double

        ' Decompose eigenvector matrix subroutine


        Eps = 0.00000001
        ITMAX = 10
        IT = 0
        NI = N * (N - 1) / 2
        KI = NI

        Do
            IT = IT + 1

            For j = 1 To N - 1
                For h = j + 1 To N
                    SCALPROD = 0.0#
                    XJLONG = 0.0#
                    XHLONG = 0.0#

                    For i = 1 To N
                        WIJ = BMATRIX(i, j)
                        WIH = BMATRIX(i, h)
                        SCALPROD = SCALPROD + WIJ * WIH
                        XJLONG = XJLONG + WIJ * WIJ
                        XHLONG = XHLONG + WIH * WIH
                    Next i

                    XLE = System.Math.Sqrt(XHLONG * XJLONG)

                    If (XLE > 0.0#) Then
                        AC = System.Math.Abs(SCALPROD / XLE)
                    Else
                        AC = 0.0#
                    End If

                    SCALPROD = SCALPROD + SCALPROD
                    WW = XJLONG - XHLONG

                    If (AC > Eps Or WW < 0.0#) Then
                        If (System.Math.Abs(WW) >= System.Math.Abs(SCALPROD)) Then
                            TETA = System.Math.Abs(SCALPROD / WW)
                            CTETA = 1.0# / System.Math.Sqrt(1.0# + TETA * TETA)
                            STETA = TETA * CTETA
                        Else
                            TETA = System.Math.Abs(WW / SCALPROD)
                            STETA = 1.0# / System.Math.Sqrt(1.0# + TETA * TETA)
                            CTETA = STETA * TETA
                        End If

                        CTETA = System.Math.Sqrt((1.0# + CTETA) / 2.0#)
                        STETA = STETA / (CTETA + CTETA)

                        If (WW < 0.0#) Then
                            Flip = CTETA
                            CTETA = STETA
                            STETA = Flip
                        End If

                        If (SCALPROD < 0.0#) Then STETA = -STETA

                        For i = 1 To N
                            Flip = BMATRIX(i, j)
                            BMATRIX(i, j) = BMATRIX(i, h) * STETA + Flip * CTETA
                            BMATRIX(i, h) = BMATRIX(i, h) * CTETA - Flip * STETA
                        Next i

                        KI = NI
                    Else
                        KI = KI - 1
                    End If
                Next h
            Next j
        Loop Until ((KI <= 0) Or (IT = ITMAX))

        For j = 1 To N
            EMATRIX(j) = 0.0#
            For i = 1 To N
                EMATRIX(j) = EMATRIX(j) + BMATRIX(i, j) * BMATRIX(i, j)
            Next i
            EMATRIX(j) = System.Math.Sqrt(EMATRIX(j))
        Next j

        For j = 1 To N
            If (EMATRIX(j) <> 0.0#) Then
                For i = 1 To N
                    BMATRIX(i, j) = BMATRIX(i, j) / EMATRIX(j)
                Next i
            End If
        Next j

    End Sub

    Sub EXPLORE(ByRef N As Integer, ByRef M As Integer, ByRef x() As Double, ByRef y() As Double, ByRef p() As Double, _
                ByRef GMATRIX() As Double, ByRef d() As Double, ByRef XINFO(,) As Double, _
                ByRef EMATRIX() As Double, ByRef CMATRIX(,) As Double, ByRef BMATRIX(,) As Double, ByRef IRESP As Integer)
        Dim j As Integer
        Dim i As Integer


        Call FILL(N, M, x, y, p, GMATRIX, d, XINFO, IRESP)

        If (IRESP) Then
            For i = 1 To N
                For j = i To N
                    CMATRIX(i, j) = XINFO(i, j) / d(i) / d(j)
                    CMATRIX(j, i) = CMATRIX(i, j)
                Next j
            Next i

            For i = 1 To N
                For j = 1 To N
                    BMATRIX(i, j) = CMATRIX(i, j)
                Next j
            Next i

            Call DECOMPOSE(N, EMATRIX, BMATRIX)
        End If
    End Sub

    Sub FILL(ByRef N As Integer, ByRef M As Integer, ByRef x() As Double, ByRef y1() As Double, ByRef p() As Double, _
            ByRef GMATRIX() As Double, ByRef d() As Double, ByRef XINFO(,) As Double, ByRef IRESP As Integer)
        Dim Delta As Double
        Dim y As Double
        Dim k As Integer
        Dim j As Integer
        Dim i As Integer

        ' Fill working matrices subroutine

        Dim DIF(M) As Double

        IRESP = 1

        For i = 1 To N
            For j = i To N
                XINFO(i, j) = 0.0#
            Next j
            GMATRIX(i) = 0.0#
        Next i

        For k = 1 To M
            Call PARDIF(N, k, x, p, DIF)
            Call USERFUNC(k, x, p, y)

            Delta = y - y1(k)

            For i = 1 To N
                For j = i To N
                    XINFO(i, j) = XINFO(i, j) + DIF(i) * DIF(j)
                Next j
                GMATRIX(i) = GMATRIX(i) + DIF(i) * Delta
            Next i
        Next k

        For i = 1 To N
            d(i) = System.Math.Sqrt(XINFO(i, i))
            If (d(i) = 0.0#) Then IRESP = 0
            For j = 1 To i - 1
                XINFO(i, j) = XINFO(j, i)
            Next j
        Next i

        System.Array.Clear(DIF, 0, DIF.Length)

    End Sub

    Sub FIND(ByRef N As Integer, ByRef EMATRIX() As Double, ByRef BMATRIX(,) As Double, ByRef d() As Double, _
            ByRef XINV(,) As Double, ByRef GMATRIX() As Double, ByRef p() As Double, ByRef P1() As Double, ByRef XMY As Double)
        Dim k As Integer
        Dim a As Double
        Dim j As Integer
        Dim i As Integer


        For i = 1 To N
            For j = i To N
                a = 0.0#
                For k = 1 To N
                    If (XMY > 0.0# Or ((XMY = 0.0#) And (EMATRIX(k) > 0.000001))) Then
                        a = a + BMATRIX(i, k) * BMATRIX(j, k) / d(i) / d(j) / (EMATRIX(k) + XMY)
                    End If
                Next k
                XINV(i, j) = a
                XINV(j, i) = a
            Next j
        Next i

        For i = 1 To N
            a = 0.0#
            For j = 1 To N
                a = a - XINV(i, j) * GMATRIX(j)
            Next j
            P1(i) = p(i) + a
        Next i
    End Sub

    Sub NLFIT(ByRef N As Integer, ByRef M As Integer, ByRef x() As Double, ByRef y() As Double, ByRef Pi() As Double, _
            ByRef PU() As Double, ByRef PL() As Double, ByRef POPT() As Double, ByRef SDev() As Double, ByRef NITER As Integer)
        Dim IRD As Integer
        Dim ICONVER As Integer
        Dim F1 As Double
        Dim INSIDE As Integer
        Dim IRESP As Integer
        Dim ISIG As Integer
        Dim ISC As Integer
        Dim FOLD As Double
        Dim f As Double
        Dim i As Integer
        Dim XMY As Double
        Dim XMYMAX As Double
        Dim Eps As Double
        Dim MAXITER As Integer

        ' Levenberg-Marquardt non-linear curve-fit subroutine

        ' Input

        '  N    = number of fitting parameters
        '  M    = number of X and Y data points
        '  X()  = array of X data points ( M rows by 1 column )
        '  Y()  = array of Y data points ( M rows by 1 column )
        '  PI() = array of initial guesses for fitting parameters
        '         ( N rows by 1 column )
        '  PU() = array of upper bounds for fitting parameters
        '         ( N rows by 1 column )
        '  PL() = array of lower bounds for fitting parameters
        '         ( N rows by 1 column )

        ' Output

        '  POPT() = array of optimized fitting parameters
        '           ( N rows by 1 column )
        '  SDEV() = array of standard deviations
        '           ( N rows by 1 column )
        '  NITER  = number of algorithm iterations


        Dim p(N) As Double
        Dim POLD(N) As Double
        Dim P1(N) As Double
        Dim EMATRIX(N) As Double
        Dim GMATRIX(N) As Double
        Dim d(N) As Double
        Dim BMATRIX(N, N) As Double
        Dim XINV(N, N) As Double
        Dim CMATRIX(N, N) As Double
        Dim XINFO(N, N) As Double

        MAXITER = MaxIterations
        Eps = 0.0001
        XMYMAX = 1000000.0#
        XMY = 0.01
        NITER = 0

        For i = 1 To N
            p(i) = Pi(i)
            POLD(i) = p(i)
        Next i

        Call SQUARESUM(M, p, x, y, f)

        FOLD = f

        Do
            ISC = 0
            Call EXPLORE(N, M, x, y, p, GMATRIX, d, XINFO, EMATRIX, CMATRIX, BMATRIX, IRESP)
            NITER = NITER + 1

            If (IRESP) Then
                While (XMY <= XMYMAX) And (ISC = 0)
                    Call FIND(N, EMATRIX, BMATRIX, d, XINV, GMATRIX, p, P1, XMY)

                    INSIDE = 1
                    For i = 1 To N
                        If (P1(i) >= PU(i) Or P1(i) <= PL(i)) Then INSIDE = 0
                    Next i
                    If (INSIDE = 1) Then
                        Call SQUARESUM(M, P1, x, y, F1)
                        'following change by LCW to prevent nonconvergences 
                        'If (F1 < f) Then
                        If F1 < f + Eps Then
                            ISC = 1
                            FOLD = f
                            f = F1
                            XMY = XMY / 100.0#
                            For i = 1 To N
                                POLD(i) = p(i)
                                p(i) = P1(i)
                            Next i
                        End If
                    End If
                    XMY = XMY * 10.0#
                End While
            End If
            Call SCHECK(N, p, POLD, f, FOLD, ISIG)
        Loop Until (NITER = MAXITER) Or (ISC = 0) Or (ISIG = 0)

        If (ISC = 1) Then Call EXPLORE(N, M, x, y, p, GMATRIX, d, XINFO, EMATRIX, CMATRIX, BMATRIX, IRESP)

        If (IRESP = 0) Then
            ICONVER = 0
        Else
            IRD = 0

            For i = 1 To N
                If (EMATRIX(i) < Eps) Then IRD = 1
            Next i

            Call FIND(N, EMATRIX, BMATRIX, d, XINV, GMATRIX, p, P1, XMY)

            If (M > N) Then
                For i = 1 To N
                    SDev(i) = System.Math.Sqrt(XINV(i, i) * f / (M - N))
                Next i
                Call SCHECK(N, P1, p, f, F1, ISIG)

                ICONVER = ISIG

            End If
        End If

        ' load optimal parameter array

        For i = 1 To N
            POPT(i) = p(i)
        Next i

        System.Array.Clear(p, 0, p.Length)
        System.Array.Clear(POLD, 0, POLD.Length)
        System.Array.Clear(P1, 0, P1.Length)
        System.Array.Clear(EMATRIX, 0, EMATRIX.Length)
        System.Array.Clear(GMATRIX, 0, GMATRIX.Length)
        System.Array.Clear(d, 0, d.Length)
        System.Array.Clear(BMATRIX, 0, BMATRIX.Length)
        System.Array.Clear(XINV, 0, XINV.Length)
        System.Array.Clear(CMATRIX, 0, CMATRIX.Length)
        System.Array.Clear(XINFO, 0, XINFO.Length)

    End Sub

    Sub PARDIF(ByRef N As Integer, ByRef k As Integer, ByRef x() As Double, ByRef p() As Double, ByRef DIF() As Double)
        Dim FDEC As Double
        Dim FINC As Double
        Dim DEL As Double
        Dim PS As Double
        Dim i As Integer
        Dim XINC As Double

        ' Numerical partial derivatives subroutine

        Dim XINCP(N) As Double
        Dim XDECP(N) As Double

        XINC = 0.01

        For i = 1 To N
            XINCP(i) = p(i)
            XDECP(i) = p(i)
        Next i

        For i = 1 To N
            PS = p(i)
            DEL = XINC * PS
            If (PS = 0.0#) Then DEL = XINC
            XINCP(i) = PS + DEL
            XDECP(i) = PS - DEL
            Call USERFUNC(k, x, XINCP, FINC)
            Call USERFUNC(k, x, XDECP, FDEC)
            DIF(i) = (FINC - FDEC) / 2.0# / DEL
            XINCP(i) = PS
            XDECP(i) = PS
        Next i

        System.Array.Clear(XINCP, 0, XINCP.Length)
        System.Array.Clear(XDECP, 0, XDECP.Length)

    End Sub

    Sub SCHECK(ByRef N As Integer, ByRef Q1() As Double, ByRef Q2() As Double, ByRef FQ1 As Double, ByRef FQ2 As Double, ByRef ISIG As Integer)
        Dim i As Integer
        Dim Eps As Double

        Eps = 0.00001

        ISIG = 0

        For i = 1 To N
            If (Q1(i) <> 0.0#) Then
                If (System.Math.Abs(1.0# - Q2(i) / Q1(i)) > Eps) Then ISIG = 1
            End If
        Next i

        If (System.Math.Abs(FQ2) < 0.0000000000000001) Then
            ISIG = 0
        ElseIf (System.Math.Abs(1.0# - FQ2 / FQ1) > Eps) Then
            ISIG = 1
        End If

    End Sub

    Sub SQUARESUM(ByRef M As Integer, ByRef p() As Double, ByRef x() As Double, ByRef y1() As Double, ByRef SQRSUM As Double)
        Dim y As Double
        Dim k As Integer

        ' Sum of squares subroutine

        SQRSUM = 0.0#

        For k = 1 To M
            Call USERFUNC(k, x, p, y)
            SQRSUM = SQRSUM + (y - y1(k)) ^ 2
        Next k
    End Sub

End Module
